      subroutine modelo2dcamadasviscoacoustico()

      use globais
      implicit none
      

      real*4 ::xmin,zmin

      real*8 :: Q
      integer :: i,j
!velocidade de propagacao no ar = 340m/s
!velocidade de propagacao na agua = 1450m/s
!massa especifica da agua = 1000




! material 1
      do j=1,int((Nz-Npml+1)/4)
	  do i=Npml,Nx-Npml+1
	      c(i,j) = 1450.
	  enddo
      enddo

! material 2
      do j=int((Nz-Npml+1)/4)+1,int(2*(Nz-Npml+1)/4)
	  do i=Npml,Nx-Npml+1
	      c(i,j) = 3000.
	  enddo
      enddo

! material 3
      do j=int(2*(Nz-Npml+1)/4)+1,int(3*(Nz-Npml+1)/4)
	  do i=Npml,Nx-Npml+1
	      c(i,j) = 4000.
	  enddo
      enddo

! material 4
      do j=int(3*(Nz-Npml+1)/4)+1,(Nz-Npml+1)
	  do i=Npml,Nx-Npml+1
	      c(i,j) = 5000.
	  enddo
      enddo


!Dentro do dominio, corrige velocidade complexa
      Q=100.
      do i=Npml,nx-Npml+1
	 do j=1,nz-Npml+1
	    c(i,j) = c(i,j) * (imag*omegas(kw))**(1./Pi * atan(1./Q))
	 enddo
      enddo

      
! prolonga as propriedades do limite do dominio para as camadas de absorcao

      do i=1,Npml-1
	Q = 1.0 + (i-1)*9.0/float(Npml-1)
	do j=1,Nz
	    if(j.le.Nz-Npml+1)then
	      c(i,j) = sqrt(c(Npml,j)*CONJG(c(Npml,j))) * (1.0+imag/(2.*Q))/sqrt(1.+1/(4.0*Q**2))
	    else
	      c(i,j) = sqrt(c(Npml,Nz-npml+1)*CONJG(c(Npml,Nz-npml+1))) * (1.0+imag/(2.*Q))/sqrt(1.+1/(4.0*Q**2))
	    endif
	enddo
      enddo


      do i=Nx - Npml+1+1 , Nx
	  Q = 1.0 + (Nx-i)*9.0/float(Npml-1)
	  do j=1,Nz
	    if(j.le.Nz-Npml+1)then
      	      c(i,j) = sqrt(c(Nx - Npml +1 ,j)*CONJG(c(Nx - Npml +1 ,j))) * (1.0+imag/(2.*Q))/sqrt(1.+1/(4.0*Q**2))
	    else
      	      c(i,j) = sqrt(c(Nx - Npml +1 ,Nz-npml+1)*CONJG(c(Nx - Npml +1 ,Nz-npml+1))) * (1.0+imag/(2.*Q))/sqrt(1.+1/(4.0*Q**2))
	    endif
	  enddo
      enddo

      do j=Nz-Npml+1+1,Nz
	  Q = 1.0 + (Nz-j)*9.0/float(Npml-1)
	  do i=1+(Nz-j),Nx-(Nz-j)
      	      c(i,j) = sqrt(c(i,Nz-Npml)*conjg(c(i,Nz-Npml)))* (1.0+imag/(2.*Q))/sqrt(1.+1/(4.0*Q**2))
	  enddo
      enddo






!  por enquanto constante
      do j=1,Nz
	  do i=1,Nx
	      rho(i,j) = 2000.
	  enddo
      enddo




      ! buoyancy e bulk modulus no grid
      do i=1,Nx
	  do j=1,Nz
	      b(i,j) = 1.0/rho(i,j)
	      K(i,j) = rho(i,j)*c(i,j)**2	
! 	      K(i,j) = rho(i,j)*c(i,j)*CONJG(c(i,j))
	  enddo
      enddo


      return
      end subroutine modelo2dcamadasviscoacoustico